function F = cal_chg_bartik(param, const, eqm)
%% Choose optimal quality after Bartik shock
        
    % Define PiB
    DB_check = ((eqm.D_H) + (param.e^param.sigma*eqm.D_F*param.shock));
    D1_check = ((eqm.D_H) + (param.e^param.sigma*eqm.D_F));
    PiB = (DB_check./D1_check).^const.gamma.*eqm.Pi1; 
        
    % Export decision
    ExportCutoffQB = (const.k0*const.lnzQ-log(param.sigma*const.gamma)+log(kron(ones(param.Nf,1),PiB)-kron(ones(param.Nf,1),eqm.Pi0))-param.mu_E)/param.sigma_E;
    ExportProbQB = normcdf(ExportCutoffQB);         
    
    % Grid search
    [qB_index, ~] = grid_search(param, const, eqm.Pi0, PiB, ExportCutoffQB, ExportProbQB); 

      
%% Calculate average wage response

    % Quality level
    q_level = const.Qgrid(eqm.q_index)'; 
    qB_level = const.Qgrid(qB_index)';
    
    % Firm wage
    wageO = const.wage_grid(eqm.q_index)';
    wageB = const.wage_grid(qB_index)';
    
%         % used in earlier estimation
%         pp_lnwage = interp1(log(const.Qgrid(eqm.nq~=0)), log(const.wage_grid(eqm.nq~=0)), 'spline', 'pp');
%         hat_lnwage = @(lnq) ppval(pp_lnwage, lnq);
%         hat_wage = @(q) exp(hat_lnwage(log(q)));   
%         wageB = hat_wage(qB_level);  
    
    
    % Average quality and wage response  
    weight = eqm.density_exp;
    avg_quality_response = sum(qB_level./q_level.*weight)/sum(weight);    
    avg_wage_response = sum(wageB./wageO.*weight)/sum(weight);
    avg_wage_response2 = sum((log(wageB)-log(wageO)).*weight)/sum(weight);

    % Average wage response in by initial quality quintile
    avg_wage_response_d = sum(kron(ones(1,5),wageB./wageO.*weight).*eqm.Qindicator,1) ...
                        ./sum(kron(ones(1,5),weight).*eqm.Qindicator,1);
    avg_wage_response_d2 = sum(kron(ones(1,5),(log(wageB)-log(wageO)).*weight).*eqm.Qindicator,1) ...
                        ./sum(kron(ones(1,5),weight).*eqm.Qindicator,1);

    % Average quality response in by initial quality quintile
    avg_quality_response_d = sum(kron(ones(1,5),qB_level./q_level.*weight).*eqm.Qindicator,1) ...
                        ./sum(kron(ones(1,5),weight).*eqm.Qindicator,1);


    
%% Changes in average wage of suppliers
  
%     avg_wgt_lnwageS_q = eqm.theta_m.*eqm.c.^(param.sigma-1)./eqm.V.*sum(const.phi_y.*const.phi_v.*repmat(eqm.Psig.*(log(param.w)+const.lnwage_grid),param.Q,1),2)'; 
%     avg_wgt_lnwageS = avg_wgt_lnwageS_q(eqm.q_index)';
%     avg_wgt_lnwageSB = avg_wgt_lnwageS_q(qB_index)';
%     avg_wageS_response = sum((avg_wgt_lnwageSB - avg_wgt_lnwageS).*weight)./sum(weight);
% 
%     avg_unwgt_lnwageS_q = 1./eqm.V.*sum(const.phi_v.*repmat(eqm.Vbar.*(log(param.w)+const.lnwage_grid),param.Q,1),2)';
%     avg_unwgt_lnwageS = avg_unwgt_lnwageS_q(eqm.q_index)';
%     avg_unwgt_lnwageSB = avg_unwgt_lnwageS_q(qB_index)';
%     avg_unwgt_wageS_response = sum((avg_unwgt_lnwageSB - avg_unwgt_lnwageS).*weight)./sum(weight);


%% Changes in average wage of customers

%     TempD = eqm.theta_v./eqm.M.*eqm.c.^(param.sigma-1).*eqm.X_m;
%     TempD(isnan(TempD))=0;
%     avg_wgt_lnwageC_q = 1./eqm.D_m.*sum(const.phi_y.*const.phi_v.*repmat((TempD.*(log(param.w)+const.lnwage_grid))',1,param.Q),1);
%     avg_wgt_lnwageC = avg_wgt_lnwageC_q(eqm.q_index)';
%     avg_wgt_lnwageCB = avg_wgt_lnwageC_q(qB_index)';
%     avg_wgt_wageC_response = sum((avg_wgt_lnwageCB - avg_wgt_lnwageC).*weight)./sum(weight);
%     
%     avg_unwgt_lnwageC_q = sum(const.phi_v.*repmat((eqm.theta_v.*(log(param.w)+const.lnwage_grid))',1,param.Q),1)./sum(const.phi_v.*repmat((eqm.theta_v)',1,param.Q),1);
%     avg_unwgt_lnwageC = avg_unwgt_lnwageC_q(eqm.q_index)';
%     avg_unwgt_lnwageCB = avg_unwgt_lnwageC_q(qB_index)';
%     avg_unwgt_wageC_response = sum((avg_unwgt_lnwageCB - avg_unwgt_lnwageC).*weight)./sum(weight);


%% Export intensity

%     rB = (eqm.D_H/param.f_vH).^(1/(param.beta_v-1))./((eqm.D_H/param.f_vH).^(1/(param.beta_v-1))+(param.e^param.sigma.*eqm.D_F.*param.shock/param.f_vF).^(1/(param.beta_v-1)));
%     DB = rB.*eqm.D_H + (1-rB).*param.e^param.sigma.*eqm.D_F.*param.shock;     
%     export_intensity = 1-rB.*eqm.D_H./DB;        


%% Return

    bartik.PiB                  = PiB;
    bartik.qB_index             = qB_index;
    bartik.avg_quality_response = avg_quality_response;
    bartik.avg_wage_response    = avg_wage_response;
    bartik.avg_wage_response_d  = avg_wage_response_d;

    %bartik.export_intensity     = export_intensity;

    %bartik.avg_wgt_wageC_response       = avg_wgt_wageC_response;
    %bartik.avg_wageS_response           = avg_wageS_response;
    %bartik.avg_unwgt_wageC_response     = avg_unwgt_wageC_response;
    %bartik.avg_unwgt_wageS_response     = avg_unwgt_wageS_response;
        
    F = bartik;

end